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Abstract. - We investigate a spinless Fermi gas trapped in a honeycomb optical lattice with 
attractive nearest-neighbor interactions. At zero temperature, mean-field theory predicts three 
quantum phase transitions, two being topological. At low interactions, the system is semi- metallic. 
Increasing the interaction further, the semi-metal destabilizes into a fully gapped superfluid. At 
larger interactions, a topological transition occurs and this superfluid phase becomes gapless, 
with Dirac-like dispersion relations. Finally, increasing again the interaction, a second topological 
transition occurs and the gapless superfluid is replaced by a different fully gapped superfluid phase. 
We analyze these different quantum phases as the temperature and the lattice filling are varied. 
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j | The discovery of high-T c superconducting materials 
[TJH] has triggered numerous studies to understand the 
possible underlying pairing mechanisms at work. The 
seminal work of Bardeen, Cooper, and Schricffer [3] has 
thus been extended to more exotic situations, collectively 
known as unconventional superconductivity, where differ- 
ent types of pairing functions, resonating valence bond 
states or topological superconductors, to cite a few, have 
been put under scrutiny gHS] - However, if progresses are 
real, the nagging question of the mechanism of high-T c su- 
perconductivity still remains unanswered to date. In this 
respect, it might be interesting to look into other related 
physical systems to get further insights or clues. Indeed, 
over the past ten years, the advances in the field of ultra- 
cold atoms loaded into optical lattices have opened the un- 
precedented opportunity to study the emergence of these 
possible exotic phases in a highly controllable and accu- 
rate way [TUHT2"] . Atomic systems are often free of many 
spurious defects plaguing condensed-matter samples which 
destroy quantum coherence. Furthermore, the interaction 
strength between atoms, relative to their tunneling ampli- 
tude, can be tuned over orders of magnitude, be it with the 
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help of Feshbach resonances [13] or by increasing the opti- 
cal lattice depth, while the lattice geometry is perfectly un- 
der command. This is exemplified by the Mott-Superfluid 
transition observed with ultracold atoms fT4"Ml6) . This 
interaction energy can even be turned attractive or re- 
pulsive, a key feature in the experimental studies of the 
celebrated BEC-BCS crossover [P7hT9] . Recently, with the 
laboratory realization in 2004 of graphene sheets [3U] , the 
honeycomb lattice has attracted a lot of attention as its 
low-energy excitations around half-filling behave as mass- 
less Weyl-Dirac fermions. This situation could be eas- 
ily mimicked with ultracold atoms [2D - I23] where differ- 
ent physical models of attractive fermions in a honeycomb 
lattice have been analyzed and different surpcrfluid states 
have been proposed [2"4T - |2"5] . In particular, different ways 
of producing nearest-neighbor interactions have been stud- 
ied in [26] for graphene and in [271I29H34] for ultracold 
gases where large values seem actually reachable for com- 
posite fermions in Fermi-Bose mixtures |35j . For instance, 
nearest-neighbor interaction strengths as large as V ~ 6i 
(where t is the nearest-neighbor hopping amplitude) have 
been reported for the 171 Yb — 174 Yb mixture at zero mag- 
netic field [35] [36]. Another promising candidate is the 
6 Li— 7 Li mixture for which V can be tuned using homonu- 
clear and heteronuclear s-wave Feshbach resonances [3"Tj . 
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In this Letter, we consider a one-component fermionic 
gas loaded on the honeycomb lattice with nearest-neighbor 
attractive interactions. Using a mean-field treatment, we 
show that this system undergoes three first-order phase 
transitions at zero temperature, two of which are topolog- 
ical, as the interaction strength is increased. In particu- 
lar the system jumps from a fully gapped supcrfluid (SF) 
phase to a gaplcss one (one with zero-energy excitation 
modes) and back again to a fully gapped SF phase. 

The honeycomb lattice consists of two shifted triangular 
sublattices, one labeled with A sites and the other with B 
sites. Each A site is connected to its three adjacent B 
sites by the vectors c a (a = 1,2,3) and the honeycomb 
diamond-shaped Bravais unit cell contains exactly one A 
site and one B site, see fig. [TJ The Fermi-Hubbard tight- 
binding Hamiltonian of our one-component system reads 




Fig. 1: The Bravais lattice of a two-dimensional regular honey- 
comb lattice is a triangular lattice with a two-point basis cell 
(grey-shaded diamond-shaped area with points A and B). The 
mean-field order parameter <5 6 C 3 is defined by the compo- 
nents S a — (bjCLi) (q = 1,2,3), i and j being nearest-neighbor 
sites connected by c a . In the paper, we set the lattice constant 
a = |c Q | to unity. 
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where t is the hopping amplitude, and V the strength 
of the nearest-neighbor interaction. The fermionic oper- 
ators a, and bj annihilate a fcrmion on A and B sites 
respectively and restricts the corresponding sums to 
nearest-neighbor sites. Throughout the paper we use t as 
the energy scale and set t = 1 hereafter. 

Dropping the Hartree-Fock terms, the mean-field 
Hamiltonian in Fourier space reads 



Hmf — 



(jk4^k + 7k&kak 
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where 7 k = Y^i e (•) denotes the grand-canonical sta- 
tistical average and k spans the honeycomb first Brillouin 



zone, a regular hexagon with area SI = 87r 2 /(3v / 3a 2 )- The 
pairing field A k reads 



V 



Ak = -Jjr X^ k ~ k ' (k-k'Sk'), 



(3) 



with N c the number of lattice unit cells. We define the 
order parameter 8 G C 3 with components S a = (bjCbi) 
(a = 1,2,3), i and j being nearest-neighbor sites con- 
nected by c„. sec fig. [T] 32,33 . Due to the [/(l) invari- 
ance of H, 6 is defined up to a global phase. As H is left 
invariant by the point group C% v , any permutation of the 
components of 8 provides another mean-field solution: the 
right and left 27r/3 rotations correspond to right and left 
cyclic permutations of the S a and the reflections about c a 
to transpositions. The permutation group S3 splits the or- 
der parameter space C 3 into a direct sum of two orthogonal 
invariant subspaces. One is spanned by Ui = -^(1,1,1) 
and corresponds to the symmetric one-dimensional irre- 
ducible representation of £3 . The other one is spanned by 
u 2 = ^=(2, —1, —1) and u 3 = ^(0, 1, —1) (this particular 
choice of basis vectors will become clear later) and corre- 
sponds to the two-dimensional irreducible representation 
of S3. Writing 8 = ^2 a r] a u a , we use the gap parameter 
norm S = \fYl a |?7a| 2 and the relative weights w a = \rj a \/S 
to quantify the strength and geometry of the SF order. We 
find 



with 



A k = -Vtf.f(k) = -Vj2vaf a , 
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We note that f*(— k) = f(k). The modulus of the gap pa- 
rameter I A k | in the Brillouin zone is depicted in figsf2Ja-c) 
for each of the geometries u Q at T = and \i = 0. We 
next diagonalize the grand-canonical mean-field Hamilto- 
nian H = Hm f — fi>N through a Bogoliubov-Valatin trans- 
formation and obtain 



n 



n c Fq + J2 *U k )4As 



k,s=± 
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The operators c ks annihilates a Bogoliubov quasi-particle 
with momentum k and index s and \i is the chemical po- 
tential (fi = at half- filling) . The excitation spectrum is 
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given by E± (k) = y X ± VY where 

X = /i 2 + | 7k | 2 + (|A k | 2 + |Al k | 2 )/2 (8) 
Y = (|A k | 2 -|Al k | 2 ) 2 /4 + 25R e (Al k A k7 2 ) 
+| 7k | 2 (V + |A k | 2 + |Al k | 2 ). 

The inversion symmetry exchanging the A and B sub- 
lattices in H leads to E s (k) = E s (— k) in T~L and to the 
innocuous change 6 — > —5, already covered by the U(l) in- 
variance. The spectrum is also invariant under A k — > Al k 
or, equivalently, under S — > S* . This symmetry reflects the 
time reversal invariance of H and implies that both S and 
its time-reversed partner 6* arc mean-field solutions of %. 
If 8 and S* cannot be matched by the U(l) symmetry, i.e. 
if S is genuinely complex, then the system exhibits spon- 
taneous time-reversal symmetry breaking superfluidity. 
The free energy per unit cell at temperature T is given 

by 
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where [3 — l/(fcsT) is the inverse temperature, fcs be- 
ing the Boltzmann constant. Fq identifies with the free 
energy per unit cell at T = 0. To compute the pairing 
field, it proves numerically convenient to directly find the 
minima of F as a function of S for given values of V, fx 
and j3 and extract the corresponding values of the order 
parameter components r) a . In our computations, the r\\ 
component was kept real but the two other complex. An 
alternative equivalent procedure is to solve the three cou- 
pled gap equations obtained from (dF/drj^iS, fx, V, (3) = 
[a = 1, 2, 3). In the thermodynamic limit, they read 
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In all our calculations, the free-energy minimum was al- 
ways found for a real 5, ruling out time- reversal symmetry 
breaking in our system. 

Wc plot in figHIa) the order parameter norm 5 at fi = 
as a function of V. At T = (blue solid line), we observe 
three discontinuous jumps signaling first-order quantum 
phase transitions. The first one occurs at V\ « 3.36, the 
second one at "V% « 7.12 and the third one at V3 « 9.15. 
We observed that, for V\ < V < V2 and up to permuta- 
tions, one of the components S a is always zero while the 
other two are always opposite of each other. For V > V2, 
and again up to permutations, two of the components 5 a 
are always identical, the remaining one being different. 
This means that, for Vi < V < V2, the order param- 
eter can always be recast in the form 8 = rj 3 113, while 
for V > V2 it always reads S = T]iUi + r]2U 2 . The three 
relative weights w a at T = and /i = are plotted in 
figCUb). We further note that the dominant geometry is 
U2 in the range V2 < V < V3 since W2 ~ 1, meaning that 



^2 = £3 ~ —8i/2. However, at V = V3, the weight of 
Ui abruptly increases from w\ ~ 0.25 to W\ ~ 0.5. For 
V > V3 , one of the S a largely dominates over the other two 
and S is essentially along c a . In figj2jd) we have plotted 
|A k | for the particular case V = 11. 

The first order nature of the transition between the 
semi-metallic and the supcrfluid phase, revealed by a jump 
in the order parameter (also for T > and /1 7^ 0) is fur- 
ther corroborated by figH] Here we have plotted the free 
energy F for interaction strength corresponding to cither 
the semi-metallic regime (continuous curves) or the super- 
fluid regime (dashed curves) against 7/3 (which is the rele- 
vant parameter to describe the transition as clear from the 
previous discussion) for different values of (3 and /1. More 
precisely the three plots correspond to (a) T = and fx = 
(b) T = and fi = 0.1 and lastly (c) f3 = 5 and fx = 0. For 
these typical cases, the free energy clearly presents two 
well separated minima for each curve: one at 773 = and 
one at rj ^ 0. The global minima determines the ground- 
state of the system while the local minima corresponds to 
a metastable state. When crossing the critical interaction 
strength (i.e. when the two minima are the same), the or- 
der parameter (773) of the groundstate will jump abruptly 
between and a finite value. Hence, for these different 
cases, including those with T > and fi ^ 0, the system 
undergoes first order phase transitions. 

A deeper insight into the various emerging solutions can 
be obtained from group theoretic arguments by expanding 
the gap equations (|10j) near S — up to second order in 
S [5]. The linear term in this Ginzburg-Landau analysis 
of the gap equations is characterized by a diagonal 3x3 
matrix Mq that is split into the previous irreducible rep- 
resentations of 53. In our case, the SF onset is obtained 
when the doubly degenerate eigenvalue of Mo is unity, thus 
selecting two dominant SF orders. At this linear level, the 
SF order can a priori develop into any superposition of 112 
and U3. However the subsequent nonlinear terms lift this 
degeneracy and only specific combinations of the two SF 
orders actually survive |S]. Assuming the SF order starts 
in only one of these dominant geometries, the question 
arises whether a sub-dominant SF order could simultane- 
ously develop, as T or V change, without the need for 
an additional phase transition. For this to happen, the 
sub-dominant order must belong to the same irreducible 
representation of the (lower) symmetry group leaving in- 
variant the excitation energies associated to the dominant 
SF order [8]. In our case, U2 is invariant under swapping 
its last two components, and so is Ui , whereas U3 flips sign. 
Hence a sub-dominant ui-SF order can only develop along 
with U2. In this scenario, the Ginzburg-Landau theory 
predicts that the coefficient 772 of the dominant geometry 
U2 vanishes like (T c — T) 1 / 2 near the critical temperature, 
whereas the coefficient 771 of the sub-dominant one Ui van- 
ishes like (T c — T) 3 / 2 [5]. We have numerically verified this 
prediction. In addition the transition between the U3-SF 
order and the mixture of the 112 and Ui SF orders must 
result from two successive second-order phase transitions 
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or from a first-order one, the latter being the one observed 
at V = V 2 . 
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Fig. 2: (color online) Plot of jAk| in the hexagonal first Bril- 
louin zone of the honeycomb lattice as obtained for (a) <5 = ui , 
(b) 5 = Xl2, (c) d = 113. The interaction strength has been 
artificially set to V = 1. The Bloch wave vector k has been 
expressed in units of k — 4n/(3a) so that the ranges in the 
hexagonal Brillouin zone are 

IW«| < V 2 and l*w/«l < 
respectively. Plot (d) is obtained from the actual minimization 
of the free-energy at V = 11 (S is a certain linear combina- 
tion of ui and 112). In this last figure the color scale has been 
modified to help visualize the pattern. 

Fig|3ja) also shows how S varies with V as T is varied. 
At j3 = 5 (red dotted line), the physics is qualitatively the 
same as at T = 0. At sufficiently high T, e.g. /? = 1 (green 
dashed line), the first and third transitions smoothen and 
become of higher order. The second transition, which im- 
plies a radical change in the geometry of 5, is always first- 
order as long as U3 persists. Above some T ((3 = 0.75), 113 
no longer appears. 

The nature of each SF order is inferred from the energy 
gap -Egap = mink (2i?_(k)), he. from the minimum en- 
ergy needed to create an excitation in the system. From 
figEt c ) j we see that at fj, = and T = the SF order 
is fully gapped below V2, gapless for V2 < V < V3 and 
again fully gapped above V3; the transitions at V% and V3 
are thus TQPTs |38l39j . The Chern number |30P2] of the 
bands is zero and, given the symmetry of the Hamiltonian, 
the zero-energy modes appear in pairs and are not topolog- 
ically protected [43]. We next verified that these TQPTs 
persist at higher T and away from half-filling. FigJ^c) 
shows £ga P at [i = for different T as V increases. For 
T as large as at f3 = 1, we still observe three distinct SF 
regions, their specific ranges in V depending on T. Fixing 




Fig. 3: (color online) (a) Plot of the order parameter strength 
8 vs V for fj, = T = (blue continuous line), j3 — 5 (red dotted 
line), P — 1.5 (green dashed line) and j3 — 0.75 (black dash- 
dotted line), (b) Plot of the geometry weights Wi (red dotted 
line), W2 (green dashed line) and W3 (blue continuous line) at 
H — T = 0. (c) Energy gap -Egap at /i = versus interaction 
strength V for T = (blue continous line), {3 = 2 (red dotted 
line), f) — 1.25 (green dashed line) and j3 = 1 (black dot-dashed 
line), (d) Excitation gap -E gap at = 2 versus interaction 
strength V for fi = (red dotted line), /x = 0.3 (green dashed 
line) and /i — 0.5 (black dot-dashed line). 



T at the experimentally attainable value j3 = 2, figEJd) 
gives -Egap for different chemical potentials. As large val- 
ues of V could be achieved [35] , the observation of these 
two TQPTs with ultracold gases seems within experimen- 
tal reach in the near future. 

FigJSJgives further evidence of the dramatic changes oc- 
curring in the excitation spectrum E_ (k) in the Brillouin 
zone when the different SF orders are achieved by increas- 
ing V at T = and fj, = 0. In the semi-metallic region, see 
figEIa), one recovers the usual honeycomb lattice tight- 
binding band structure when V — > 0. FigGJb) gives the 
spectrum for the fully gapped SF phase at V = 5. At 
V = 8, a gapless SF order develops and -E_(k) presents 
two zeros with linear dispersion, see figlSJc). The tips of 
these cones move closer to each other when V is increased 
but do not merge at T = 0. Instead, the spectrum changes 
abruptly at V3 and displays a single nonzero minimum, sig- 
naling again a fully gapped SF order. This is shown in fig[5] 
for V = 11. At higher T, this transition becomes smooth. 
The two cones merge into one single zero energy minimum 
with linear dispersion in one direction and quadratic in the 
orthogonal one. A similar-looking topological transition 
is obtained for non-interacting fermions in graphene when 
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Fig. 4: (color online) Plot of the free energies F versus 773 
(771 = 112 = 0) below the transition at V\ (semi- metallic regime) 
and above Vi (superfluid regime), (a) T = 0, fj, = 0, V = 3.3 
(continuous line) and V = 3.4 (dashed line); (b) T = 0, fi = 
0.1, V = 3.26 (continuous line) and V = 3.38 (dashed line); 
(c) j3 = 5, fi = 0, V = 3.35 (continuous line) and V = 3.42 
(dashed line). The local minimum at r\ 7^ (continuous line) 
becomes a global one (dashed line) leading to a jump of the 
order parameter. Therefore the transition at Vi is of the first- 
order also for a moderate increase of the temperature as well 
as a departure from half-filling. 



the nearest-neighbor hopping amplitudes are imbalanccd 
[HI [45]. The energy gap -E gap then increases smoothly 
from zero when V is further increased and the system en- 
ters again a fully gapped SF order. 

As a conclusion, we have observed and studied two 
topological quantum phase transitions in a one-component 
Fermi gas loaded in a honeycomb lattice with nearest- 
neighbour attractive interactions. The corresponding SF 
features are robust against moderate changes in the tem- 
perature and in the chemical potential. These topological 
transitions should be within the reach of ultracold gases 
experiments. The different pairing geometries achieved by 
the system could be analyzed by observing the momentum 
distribution of the expanding gas after release from the 
trap. 

To estimate the critical temperature of the transi- 
tion (Kosterliz-Thouless-like) between the gas of paired 
fermions and the superfluid phase, one must take into ac- 
count the gaussian fluctuations beyond the saddle-point 
mean-field approximation 24,38,39 . This will be under- 
taken in a future work. 



E 




k fK 



in 
J 
-1/2 -1/2 1/2 



.i /2 -1/2 k IK 

k Ik * 

y 




Fig. 5: (color online) Plots of E-(k) in the Brillouin zone 
at /j, = T = 0. (a) Honeycomb lattice tight-binding band 
structure |-yic| at V = 0. (b) Fully-gapped SF phase at V — 5. 
(c) Gapless SF phase at V = 8. (d) Fully gapped SF phase at 
V = 11. The Bloch wave vector k has been expressed in units 
of k = 4-7r/(3a) so that the ranges in the hexagonal Brillouin 
zone are |fca,//t| < 1/2 and \k y /K\ < l/s/3 respectively. 
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